function [phi,PMarkov] = dreumProbWTP2(H,K,Np,quota,sigma)

% Calculate success probabilities for each lottery (PHI), and Markov transition
% matrix (PMarkov)


% Calculate probability of success for each hunt---------------------------
phi = zeros(H,K); % Initialize array for storing probability of success


for idx1 = 2:H
        
    for idx2 = 1:K
        
        % Calculate applicants with k or more preference points that apply 
        % to hunt h
        applicantsp = Np(idx2:K)*sigma(idx1,idx2:K)';
        
        % Calculate applicants with k + 1 or more preference points that
        % apply to hunt h
        if idx2 < K
            
            applicantspPlus1 = Np(idx2+1:K)*sigma(idx1,idx2+1:K)';
            
        else
            
            applicantspPlus1 = Np(idx2:K)*sigma(idx1,idx2:K)';
            
        end
        
        % Calculate first-round success probabilities
        if applicantsp <= quota(idx1-1)
            
            phi(idx1,idx2) = 1; % Win with certainty if more 
                                  % permits than applicants
                                  
        elseif ((applicantsp > quota(idx1-1)) && ...
                    (applicantspPlus1 <= quota(idx1-1))) 
                           
            % Simple lottery for permits that remain after higher permit
            % holders get theirs
            phi(idx1,idx2) = (quota(idx1-1) - applicantspPlus1)...
                    /(Np(idx2)*sigma(idx1,idx2));
                
        elseif (applicantspPlus1 >= quota(idx1-1) && idx2 == K)
            
            phi(idx1,idx2) = quota(idx1-1)/(Np(idx2)*sigma(idx1,idx2));
                    
        end
            
    end
    
end
    
phi = [phi; ones(1,K)];


% Calculate Markov transition matrix
PMarkov = zeros(K,K,H+1);
% PMarkov = cell(1,2);

for idx1 = 1:H+1 % idx1 = 1 == pp-only option

    if idx1 < H+1
    
        for idx2 = 1:K
        
            % Assign success probabilities
            P = phi(idx1,idx2);
                    
            if idx2 < K
            
                if idx1 == 1 % So that we are on hunt 0 (pp-only option)
                
                    % Automatically win another point if pp-only option chosen
                    PMarkov(idx2,idx2+1,idx1) = 1;
            
                else
                
                    % Transition to zero pts if successfully drawn
                    PMarkov(idx2,1,idx1) = P;
                
                    % Transition to k + 1 pts if not drawn
                    PMarkov(idx2,idx2+1,idx1) = 1 - P;
                
                end
            
%             elseif idx3 == 1
%                 
%                 % Transition to zero pts if successfully drawn
%                 x(K(idx3),1,idx1,idx3) = P;
%                 
%                 % Transition to K + 1 pts if not drawn (yr 1 only)
%                 x(K(idx3),K(idx3)+1,idx1,idx3) = 1 - P;
                
            else
            
                % Transition to zero pts if successfully drawn
                PMarkov(K,1,idx1) = P;
                
                % Stay at K pts if not drawn
                PMarkov(K,K,idx1) = 1 - P;
                
            end
            
        end
        
    else
        
        PMarkov(:,:,idx1) = eye(K);
    end
        
end
    
    
% % Calculate shares of applicants for each site choice (note:
% % there are 22 first-round choices, plus the preference point only option => 
% % 23 site choice combinations)
% sigma{1} = numPatch.*ones(H,K(1));
% sigma{2} = numPatch.*ones(H,K(2)); % If no one with k pts applied to site h, 
%                                      % set share to numPatch do avoid undefined 
%                                      % variable (since ln(0) = undef)
% 
% for idx3 = 1:2
%                                      
%     for idx1 = 1:H
%     
%         for idx2 = 1:K(idx3)
%         
%             % Identify applicants who applied to hunt h w/ k preference points
%             huntCond = (data(:,3,idx3) == idx1) & (data(:,4,idx3) == idx2 - 1);
%         
%             % Identify applicants with k pts
%             pointCond = (data(:,4,idx3) == idx2 - 1) & (data(:,3,idx3) ~= 24); 
%             
%             % Calculate share of hunters who apply to each sige choice combo
%             sigma{idx3}(idx1,idx2) = size(data(huntCond,:,idx3),1)/size(data(pointCond,:,idx3),1);
%         
%         end
%     
%     end
%     
%     sigma{idx3}(sigma{idx3}==0) = numPatch;
%     
% end   
